Setup

library(sp)
library(raster)
library(R.matlab)
library(plyr)
library(data.table)
library(maptools)
library(rgdal)
library(spatstat) 

Data Proccessing

setwd("~/Desktop/Fall_2018_Classes/GIS_713/Final_Project")
CHIPS_df <- read.table("chips.csv", header = TRUE, row.names=NULL, sep=",")
CHIPS_df <- CHIPS_df[!CHIPS_df$lat < 500000.00, ]
CHIPS_df <- CHIPS_df[!CHIPS_df$long < 50000.00, ]
# Correcting odd lat coords
CHIPS_df$lat <- CHIPS_df$lat / 10000
CHIPS_df$long <- CHIPS_df$long / 10000
# And a function to shift vectors conviniently:
shift.vec <- function (vec, shift) {
  if(length(vec) <= abs(shift)) {
    rep(NA ,length(vec))
  }else{
    if (shift >= 0) {
      c(rep(NA, shift), vec[1:(length(vec)-shift)]) }
    else {
      c(vec[(abs(shift)+1):length(vec)], rep(NA, abs(shift))) } }
  }
# Calculating distances between successive positions and the respective speed in this segment.
# Shifting vectors for lat and lon so that each row also contains the next position.
CHIPS_df$lat.p1 <- shift.vec(CHIPS_df$lat, -1)
CHIPS_df$lon.p1 <- shift.vec(CHIPS_df$long, -1)
# Calculating distances (in metres) using the function pointDistance from the ‘raster’ package.
CHIPS_df$dist.to.prev <- apply(CHIPS_df, 1, FUN = function (row) {
  pointDistance(c(as.numeric(row["lat.p1"]), as.numeric(row["long.p1"])), 
                c(as.numeric(row["lat"]), as.numeric(row["long"])), 
                lonlat = T)
})
# Transforming the column ‘time’ so that R knows how to interpret it.
CHIPS_df$time_new <- strptime(CHIPS_df$initial_time_stamp_mat,
                              format="%m/%d/%Y %H:%M")
# Shift the time vector, too.
CHIPS_df$time.p1 <- shift.vec(CHIPS_df$time_new, -1)
# Calculating number of seconds between two positions.
CHIPS_df$time.diff.to.prev <- as.numeric(difftime(CHIPS_df$time.p1, 
                                                  CHIPS_df$time_new))
# Calculating metres per seconds, kilometres per hour, and two LOWESS smoothers to get rid of some noise.
CHIPS_df$speed.m.per.sec <- CHIPS_df$dist.to.prev / CHIPS_df$time.diff.to.prev
CHIPS_df$speed.km.per.h <- CHIPS_df$speed.m.per.sec * 3.6
CHIPS_df$speed.km.per.h <- ifelse(is.na(CHIPS_df$speed.km.per.h), 0, 
                                  CHIPS_df$speed.km.per.h)
CHIPS_df$lowess.speed <- lowess(CHIPS_df$speed2, f = 0.2)$y
CHIPS_df$lowess.alt <- lowess(CHIPS_df$altitude, f = 0.2)$y
CHIPS_df$lowess.conduct <- lowess(CHIPS_df$conductance_z, f = 0.2)$y

setnames(CHIPS_df, "long", "lon")
pt1 <- CHIPS_df[CHIPS_df$participant == 1, ]
pt2 <- CHIPS_df[CHIPS_df$participant == 2, ]
pt3 <- CHIPS_df[CHIPS_df$participant == 3, ]
pt4 <- CHIPS_df[CHIPS_df$participant == 4, ]
pt5 <- CHIPS_df[CHIPS_df$participant == 5, ]
pt6 <- CHIPS_df[CHIPS_df$participant == 6, ]
pt7 <- CHIPS_df[CHIPS_df$participant == 7, ]
pt8 <- CHIPS_df[CHIPS_df$participant == 8, ]
pt9 <- CHIPS_df[CHIPS_df$participant == 9, ]
pt10 <- CHIPS_df[CHIPS_df$participant == 10, ]
pt11 <- CHIPS_df[CHIPS_df$participant == 11, ]

Exploratory Data Visualizations

GPS

# Now, let’s plot all the stuff!
# Plot elevations and smoother
layout(matrix(1:3, nrow=3))
plot(CHIPS_df$altitude, type = "l", bty = "n", xaxt = "n", lwd= 3, 
     ylab = "Elevation", 
     xlab = "", col = "grey60")
lines(CHIPS_df$lowess.alt, col = "green", lwd = 3)
abline(h = mean(CHIPS_df$altitude), lty = 2, lwd = 3, col = "green")
legend(x="bottomright", legend = c("GPS elevation", "LOWESS elevation", 
                                "Mean elevation"),
       col = c("grey60", "green", "green"), lwd = c(2,4,2), lty = c(1,2,2),
       bty = "n")
# Plot speeds and smoother
plot(CHIPS_df$speed2, type = "l", bty = "n", lwd= 3, xaxt = "n", 
     ylab = "Speed (km/h)", xlab = "", col = "grey60")
lines(CHIPS_df$lowess.speed, col = "red", lwd = 3)
abline(h = mean(CHIPS_df$speed2), lty = 2, lwd = 3, col = "red")
legend(x="topright", legend = c("GPS speed", "LOWESS speed", 
                                   "Mean speed"),
       col = c("grey60", "red", "red"), lwd = c(2,4,2), lty = c(1,2,2), 
       bty = "n")
# Plot conductnace and smoother
plot(CHIPS_df$conductance_z, type = "l", bty = "n", lwd= 3, xaxt = "n", 
     ylab = "Skin Conductance", xlab = "", col = "grey60")
lines(CHIPS_df$lowess.conduct, col = "blue", lwd = 3)
abline(h = mean(CHIPS_df$conductance_z), lty = 2, lwd = 3,  col = "blue")
legend(x="topright",
       legend = c("Conductance", "LOWESS conductance", "Mean conductance"),
       col = c("grey60", "blue", "blue"), lwd = c(2,4,2), lty = c(1,2,2),
       bty = "n")

par(mar=c(5, 4, 4, 2) + 0.1)

Skin Conductance

Plotting the elevation and timestamp of each waypoint using ggplot, allows the visualization of the cycling route between Tilburg and Waalwijk. As a preparatory step, the ymd_hms() function from the lubridate library is used to convert the string representating the timestamp to a proper R time-object. To not confuse ggplot, the SpatialPointsDataFrame-object is not passed directly, but converted to a regular dataframe with as.data.frame():

if(!require(lubridate)) install_github("rstudio/lubridate")
if(!require(ggplot2)) install_github("rstudio/ggplot2")
if(!require(gridExtra)) install_github("rstudio/gridExtra")
# plot of time and elevation, colored by skin conductance
time_ele_conduct_plot <- ggplot(as.data.frame(CHIPS_df), # convert to regular dataframe
            aes(x=time, y=altitude, color = conductance_z)) +
            scale_color_gradient(low="blue", high="red") +
            geom_point(alpha = 0.01, size = 2) + 
            labs(x='Cycling time', y='Elevations (meters)')
# plot of time and speed, colored by skin conductance
time_speed_conduct_plot <- ggplot(as.data.frame(CHIPS_df), # convert to regular dataframe
            aes(x=time, y=speed2, color = conductance_z)) +
            scale_color_gradient(low="blue", high="red") +
            geom_point(alpha = 0.08, size = 2) + 
  labs(x='Cycling time', y='Speed (km/h)')
grid.arrange(time_ele_conduct_plot, time_speed_conduct_plot, nrow=2)

Spatial Data Processing

shp_dsn <- "~/Desktop/Fall_2018_Classes/GIS_713/Final_Project/NL006L3_TILBURG/Shapefiles"
landcover <- readOGR(path.expand(shp_dsn), 'NL006L3_TILBURG_UA2012')
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/garrettmillar/Desktop/Fall_2018_Classes/GIS_713/Final_Project/NL006L3_TILBURG/Shapefiles", layer: "NL006L3_TILBURG_UA2012"
## with 12086 features
## It has 11 fields
# Projection
landcover <- spTransform(landcover, CRS("+proj=longlat +datum=WGS84"))
# Conversion into SpatialPoints
coordinates(CHIPS_df) <- ~ lon + lat
coordinates(pt1) <- ~ lon + lat
coordinates(pt2) <- ~ lon + lat
coordinates(pt3) <- ~ lon + lat
coordinates(pt4) <- ~ lon + lat
coordinates(pt5) <- ~ lon + lat
coordinates(pt6) <- ~ lon + lat
coordinates(pt7) <- ~ lon + lat
coordinates(pt8) <- ~ lon + lat
coordinates(pt9) <- ~ lon + lat
coordinates(pt10) <- ~ lon + lat
coordinates(pt11) <- ~ lon + lat
# Setting default projection
proj4string(CHIPS_df) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt1) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt2) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt3) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt4) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt5) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt6) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt7) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt8) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt9) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt10) <- CRS('+proj=longlat +datum=WGS84')
proj4string(pt11) <- CRS('+proj=longlat +datum=WGS84')

Web-mapping

library(leaflet)
require(pals)
conduct.pal <- colorNumeric (c("dodgerblue4", "slategray2", "red3"), 
                             pt1$conductance_z)
m <- leaflet() %>%
  # Add tiles
  addProviderTiles("Esri.WorldTopoMap", group = "Topographical") %>%
  addProviderTiles("OpenStreetMap.Mapnik", group = "Road map") %>%
  addProviderTiles("Esri.WorldImagery", group = "Satellite") %>%
  addCircles (data=pt1, group='Participant 1', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt1$conductance_z),
              opacity = 0.2, color = conduct.pal(pt1$conductance_z)) %>%
  addCircles (data=pt2, group='Participant 2', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt2$conductance_z),
              opacity = 0.2, color = conduct.pal(pt2$conductance_z)) %>%
  addCircles (data=pt3, group='Participant 3', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt3$conductance_z),
              opacity = 0.2, color = conduct.pal(pt3$conductance_z)) %>%
  addCircles (data=pt4, group='Participant 4', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt4$conductance_z),
              opacity = 0.2, color = conduct.pal(pt4$conductance_z)) %>%
  addCircles (data=pt5, group='Participant 5', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt5$conductance_z),
              opacity = 0.2, color = conduct.pal(pt5$conductance_z)) %>%
  addCircles (data=pt6, group='Participant 6', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt6$conductance_z),
              opacity = 0.2, color = conduct.pal(pt6$conductance_z)) %>%
  addCircles (data=pt7, group='Participant 7', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt7$conductance_z),
              opacity = 0.2, color = conduct.pal(pt7$conductance_z)) %>%
  addCircles (data=pt8, group='Participant 8', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt8$conductance_z),
              opacity = 0.2, color = conduct.pal(pt8$conductance_z)) %>%
  addCircles (data=pt9, group='Participant 9', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt9$conductance_z),
              opacity = 0.2, color = conduct.pal(pt9$conductance_z)) %>%
  addCircles (data=pt10, group='Participant 10', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt10$conductance_z),
              opacity = 0.2, color = conduct.pal(pt10$conductance_z)) %>%
  addCircles (data=pt11, group='Participant 11', stroke = T, radius = 80, 
              fillOpacity = 0.2, fillColor = conduct.pal(pt11$conductance_z),
              opacity = 0.2, color = conduct.pal(pt11$conductance_z)) %>%
  # Layers control
  addLayersControl(position = 'bottomright',
                   baseGroups = c("Topographical", "Road map", "Satellite"),
                   overlayGroups = c("Participant 1", "Participant 2", 
                                     "Participant 3", "Participant 4", 
                                     "Participant 5", "Participant 6",
                                     "Participant 7", "Participant 8", 
                                     "Participant 9", "Participant 10", 
                                     "Participant 11"),
                   options = layersControlOptions(collapsed = FALSE)) %>%
  hideGroup(c("Participant 2", "Participant 3", "Participant 4", "Participant 5", 
              "Participant 6", "Participant 7", "Participant 8", "Participant 9",
              "Participant 10", "Participant 11")) %>%
  addLegend(values = pt1$conductance_z, pal = conduct.pal, 
            opacity = 1, title = "Skin Conductivity", position = "bottomleft")
m

Study Area

urban <-  c("Discontinuous dense urban fabric (S.L. : 50% -  80%)",  
            "Industrial, commercial, public, military and private units",
            "Discontinuous low density urban fabric (S.L. : 10% - 30%)", 
            "Isolated structures",
            "Continuous urban fabric (S.L. : > 80%)",
            "Discontinuous very low density urban fabric (S.L. : < 10%)",               
            "Discontinuous medium density urban fabric (S.L. : 30% - 50%)",  
            "Construction sites",  
            "Mineral extraction and dump sites")
natural <- c("Arable land (annual crops)",  
             "Pastures",   
             "Forests",    
             "Land without current use",
             "Green urban areas", 
             "Permanent crops (vineyards, fruit trees, olive groves)", 
             "Sports and leisure facilities",   
             "Herbaceous vegetation associations (natural grassland, moors...)",
             "Open spaces with little or no vegetation (beaches, dunes, bare rocks, glaciers)",
             "Wetlands")
roads  <- c("Other roads and associated land",
            "Railways and associated land",
            "Fast transit roads and associated land")
name = landcover$ITEM2012
landcover$group <- with(landcover, ifelse(name %in% urban, "urban",
                                          ifelse(name %in% natural, 
                                                 "natural", "roads")))
col = terrain.colors(3)
plot(landcover, col = col, col.regions = landcover$group,
                 edge.col = "transparent", axes = F,
                 colorkey = list(space = "bottom", height = 0.5, 
                                 width = 0.7),
                 main = "Study Area", sub = "Types of Landcover",
                par.settings = list(axis.line = list(col = 'transparent')))
legend("bottomleft", title = "Landcover", fill = col, 
       legend = c("Urban", "Natural", "Roads"))

Spatial Analysis

Mapping Cycling Routes

# clipping landcover polygon to cycling route
neth_clipped <- as(crop(landcover, CHIPS_df), "SpatialPolygonsDataFrame")
col = terrain.colors(3)
conduct.pal <- colorNumeric (c("dodgerblue4", "slategray2", "red3"), 
                             pt1$conductance_z)
layout(matrix(1:6, nrow=1))
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt1, col = conduct.pal(pt1$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt2, col = conduct.pal(pt2$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt3, col = conduct.pal(pt3$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt4, col = conduct.pal(pt4$conductance_z), cex = 1.5, pch = 20 )
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt5, col = conduct.pal(pt5$conductance_z), cex = 1.5, pch = 20)
plot(neth_clipped, col = col, col.regions = neth_clipped$group,
     edge.col = "transparent", axes = F, 
     colorkey = list(space = "bottom", height = 0.5, width = 0.7),
     main = "", 
     par.settings = list(axis.line = list(col = 'transparent')))
points(pt6, col = conduct.pal(pt6$conductance_z), cex = 1.5, pch = 20)

Skin Conductance by Landcover Types

require(spatialEco)
require(sp)
neth.pts.poly <- point.in.poly(CHIPS_df, neth_clipped)
# Number of points in each landcover group
Conduct_Points_per_group <- tapply(neth.pts.poly@data$conductance_z, 
                                  neth.pts.poly@data$group, 
                                  FUN=length)
# Mean conductance in each landcover group
Conduct_Mean_per_group <- tapply(neth.pts.poly@data$conductance_z, neth.pts.poly@data$group, FUN=mean)
Conduct_Mean_per_group <- round(Conduct_Mean_per_group, 2)
conductance_groups <- data.frame(t(rbind("Skin Conductance Points (N)"=Conduct_Points_per_group,
                                "Mean Skin Conductance"=Conduct_Mean_per_group)))
# Number of points in each polygon
Conduct_Points_per_Poly <- tapply(neth.pts.poly@data$conductance_z, neth.pts.poly@data$ITEM2012, 
                                  FUN=length)
# Mean conductance in each polygon
Conduct_Mean_per_Poly <- tapply(neth.pts.poly@data$conductance_z, neth.pts.poly@data$ITEM2012, FUN=mean)
Conduct_Mean_per_Poly <- round(Conduct_Mean_per_Poly, 2)
conductance_polys <- data.frame(t(rbind("Skin Conductance Points (N)"=Conduct_Points_per_Poly,
                                "Mean Skin Conductance"=Conduct_Mean_per_Poly)))
library(tableHTML)
# Table: Number of points and means (landcover groups)
conductance_groups[is.na(conductance_groups)] <- 0
conductance_groups_table <- conductance_groups %>%
  tableHTML( border = 2,
          rownames = TRUE, 
          headers = c("Sampled Points (N)", "Skin Conductance (M)"), 
          second_headers = list(c(2, 2), c('Landcover Class', 'Statistic'))) %>%
  add_css_second_header(css = list(c('background-color', 'color', 'height'),
                                   c('#C0C0C0', 'black', '50px')),
                        second_headers = 1:2) %>%
    add_css_row(css = list('background-color', '#f2f2f2'), rows = odd(3:6)) %>%
  add_css_conditional_column(conditional = ">",
                             value = 4000,
                             css = list('background-color', "lightcoral"),
                             columns = c("Mean Skin Conductance")) %>%
  add_css_conditional_column(conditional = "<=",
                             value = -1000,
                             css = list('background-color', "lightsteelblue"),
                             columns = c("Skin Conductance (M)"))
# Table: Number of points and means (all landcover types)
conductance_polys[is.na(conductance_polys)] <- 0
conductance_polys_table <- conductance_polys %>%
  tableHTML( border = 2,
          rownames = TRUE, 
          headers = c("Sampled Points (N)", "Skin Conductance (M)"), 
          second_headers = list(c(1, 3), c('Landcover Class', 'Statistic'))) %>%
  add_css_second_header(css = list(c('background-color', 'color', 'height'),
                                   c('#C0C0C0', 'black', '50px')),
                        second_headers = 1:2) %>%
    add_css_row(css = list('background-color', '#f2f2f2'), rows = odd(3:25)) %>%
  add_css_conditional_column(conditional = ">",
                             value = 4000,
                             css = list('background-color', "lightcoral"),
                             columns = c("Skin Conductance (M)")) %>%
  add_css_conditional_column(conditional = "<=",
                             value = -1000,
                             css = list('background-color', "lightsteelblue"),
                             columns = c("Skin Conductance (M)"))

conductance_groups_table
Landcover Class Statistic
Sampled Points (N) Skin Conductance (M)
natural 33024 3250.58
roads 46954 648.91
urban 28333 1235.15
conductance_polys_table
Landcover Class Statistic
Sampled Points (N) Skin Conductance (M)
Arable land (annual crops) 6480 3134.32
Construction sites 0 0
Continuous urban fabric (S.L. : > 80%) 3995 1574.56
Discontinuous dense urban fabric (S.L. : 50% - 80%) 4004 -2737.85
Discontinuous low density urban fabric (S.L. : 10% - 30%) 2764 5136.59
Discontinuous medium density urban fabric (S.L. : 30% - 50%) 1600 -2537.95
Discontinuous very low density urban fabric (S.L. : < 10%) 40 -0.3
Fast transit roads and associated land 0 0
Forests 9408 2292.75
Green urban areas 268 -1840.75
Herbaceous vegetation associations (natural grassland, moors...) 2668 540.78
Industrial, commercial, public, military and private units 15522 1851.4
Isolated structures 408 1944.19
Land without current use 1224 1316.95
Mineral extraction and dump sites 0 0
Open spaces with little or no vegetation (beaches, dunes, bare rocks, glaciers) 0 0
Other roads and associated land 43466 421.57
Pastures 5332 867.85
Permanent crops (vineyards, fruit trees, olive groves) 340 10679.79
Railways and associated land 3380 3439.65
Sports and leisure facilities 7304 7481.75
Water 108 4802.31
Wetlands 0 0

Space-time Cube

library(OpenStreetMap)
map <- openmap(as.numeric(c(max(pt1$lat), min(pt1$lon))),
               as.numeric(c(min(pt1$lat), max(pt1$lon))), 
               type = "stamen-terrain")
transmap <- openproj(map, projection = proj4string(pt1))
map3d <- function(map, ...){
  if(length(map$tiles)!=1){stop("multiple tiles not implemented") }
  nx = map$tiles[[1]]$xres
  ny = map$tiles[[1]]$yres
  xmin = map$tiles[[1]]$bbox$p1[1]
  xmax = map$tiles[[1]]$bbox$p2[1]
  ymin = map$tiles[[1]]$bbox$p1[2]
  ymax = map$tiles[[1]]$bbox$p2[2]
  xc = seq(xmin,xmax,len=ny)
  yc = seq(ymin,ymax,len=nx)
  colours = matrix(map$tiles[[1]]$colorData,ny,nx)
  m = matrix(0,ny,nx)
  surface3d(xc,yc,m,col=colours, ...) }
library(RColorBrewer)
bp = brewer.pal(11,"RdBu")
library(colourschemes)
cs = rampInterpolate(pt1$conductance_z, rev(bp))
pt1$conduct_pal <- cs(pt1$conductance_z)
plot3d(pt1$lon, pt1$lat, pt1$time, xlab="Longitude", 
       ylab="Latitude", zlab="Time", type = "s", 
       col = pt1$conduct_pal, size = 2.5) 
map3d(transmap)

You must enable Javascript to view this page properly.